// Локальные
#include "linal/bcg.h"
#include "linal/gauss.h"

namespace linal
{
    //============================================================================================================================
    //============================================================================================================================
    double BCG::solve(const CSRMatrix &matrix, const Vector &f, Vector &x)
    {
        double alpha, beta;	//Вспомогательные константы метода.
        double norm_f;		//Норма вектора правой части (Чтобы не пересчитывать много раз).

        // Нахождение транспонированной матрицы.
        CSRMatrix transposed_matrix = matrix.transposed();

        // Вычисление нормы f
        norm_f = f.normSquared();

        // Нахождение r нуле-вого.
        m_R = matrix*x;
        m_R = f - m_R;

        // Нахождение z нулевого.
        m_P = m_R;

        // Нахождение p нулевого.
        m_Z = m_R;

        // Нахождение s нулевого.
        m_S = m_R;

        // Итерационный шаг.
        double residual = std::numeric_limits<double>::max();
        for(int cur_iter = 0; cur_iter < m_MaxNumOfIter && residual > m_SufficientResidual*m_SufficientResidual; cur_iter++)
        {
            m_Az = matrix*m_Z;
            m_ATs = transposed_matrix*m_S;

            double pr_scalar = m_P*m_R;

            // Нахождение alpha i-го.
            alpha = pr_scalar/(m_S*m_Az);
			if(alpha < 1e-10)
				break;

            // Нахождение x i-го.
            x += alpha*m_Z;

            // Нахождение r i-го.
            m_R -= alpha*m_Az;

            // Нахождение p i-го.
            m_P -= alpha*m_ATs;

            // Нахождение beta i-го.
            beta = (m_P*m_R)/(pr_scalar);

            // Нахождение z i-го.
            m_Z = m_R + beta*m_Z;

            // Нахождение s i-го.
            m_S = m_P + beta*m_S;

            // Нахождение невязки.
            residual = m_R.normSquared()/norm_f;
        }
		return qSqrt(residual);
    }
}
